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Vanadium dioxide (VO2) undergoes a first or- 
der metal insulator transition (MIT) at 340 K 
[1]. At high temperature, the crystal structure 
is metallic with the rutile structure (R), while 
it transforms to the monoclinic (Mi) phase and 
becomes insulating below the transition temper- 
ature. Experimental evidence [2] suggests that 
the insulator is obtained due to strong correla- 
tions in this material, however density functional 
theory (DFT) predicts the Mi phase to be metal- 
lic [3, 4]. In this work, we develop and carry out 
state of the art linear scaling DFT calculations re- 
fined with non-local dynamical mean-field theory 
(DMFT). We validate our approach by applying 
this cutting edge methodology to the Mi phase 
of VO2. We find that a Peierls-Mott instability 
is responsible for the insulating Mi phase, which 
furthermore survives a large degree of disorder 
accounting for the MIT observed [5] when no 
long-range crystallographic order is present. The 
Peierls-Mott instability involves two electrons per 
vanadium atom in ti g orbitals and reconciles the 
Peierls picture with recent soft x-ray absorption 
spectroscopy (XAS) [6], which points towards a 
breaking of the one electron per 3d orbital pic- 
ture suggested early by Goodenough [7]. Finally, 
since fluctuations in VO2 thin film deposition may 
cause departure from the ideal stoichiometry [8], 
we discuss how the MIT is affected by oxygen 
vacancies. 

The nature of the metal-insulator (MIT) transition in 
VO2 has been long discussed, with particular emphasis 
placed on the role of electron correlations in forming the 
charge gap. Photoemission experiments give strong ev- 
idence for strong electron-electron and clcctron-phonon 
coupling in VO2 [2] , suggesting that this compound is an 
archetypal Mott insulator. An alternative point of view is 
that the low-temperature phase of VO2 may constitute 
a band (Peierls) insulator, where the crystal distortion 
with the V-V dimerization splits the ai g bonding band. 
Lastly one should consider a charge transfer insulator ex- 
hibiting a strong mass renormalisation. The purpose of 



this paper is to disentangle these competing pictures. 

The Peierls picture was supported from a first- 
principles perspective using LDA+GW calculations, 
where the authors found that dynamical effects com- 
bined with off-diagonal matrix elements in the self-energy 
opened a gap [4] , although the value of the fundamental 
gap was there predicted to be almost zero and thus well 
below the experimental value of 0.6 eV [9]. Very recently, 
a model Hamiltonian approach using cluster DMFT ap- 
plied to a three band Hamiltonian for the t2 9 orbitals has 
been shown to successfully capture the insulating nature 
of the Mi phase [10, 11], and the authors found a charge 
gap of 0.6 eV, in very good agreement with experiment 
[9]. Hence, V0 2 is, in the latter view, not a conventional 
Mott insulator. Instead, the formation of dynamical V-V 
singlet pairs due to strong Coulomb correlations is nec- 
essary to trigger the opening of a Peierls gap. We note, 
however, that in Ref. [10] the vanadium 3d subshell is 
occupied by a single electron (0.8 electrons for the &\ g 
with only 0.1 electrons remaining in each of the e ng or- 
bitals). A general problem with model Hamiltonian ap- 
proaches, recently pointed out in Ref. [12], is that the 
3d orbital density is very much affected by the orbital 
subset projection used in the calculations. In particular, 
it has been shown recently using XAS measurements [6] 
that the states of V0 2 are not well characterized by a 
single dominant ionic configuration, rather exhibiting a 
distributed orbital character, suggesting room for correc- 
tion of Goodenough's ionic picture of VO2. 

The key issues that we address in this work are: (1) Is 
the 3d 1 ionic picture of Goodenough valid and how many 
electrons are involved in the orbital selection process; (2) 
Can Mott correlations alone drive V0 2 to an insulator, 
and what is the minimal local repulsion Ud necessary 
to localize the charge, i.e., the Zaanen-Sawatzky- Allen 
(ZSA) boundary [13, 14]; (3) How is the ZSA boundary 
affected by other localization processes, such as the An- 
derson charge localization induced by disorder, and can 
we find an insulator for a combination of realistic disor- 
der and Coulomb repulsion; (4) Are non-local corrections 
to the self energy (the Peierls mechanism) an essential in- 
gredient to trigger the gap opening for reasonable local 
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repulsion Ug, and is the latter insulating phase stable 
against external perturbations such as disorder. 

To address these points, we move beyond the model 
Hamiltonian approach and investigate the effect of corre- 
lations in a disordered prototype for the metal-insulator 
transition (MIT) in VO2 from the ab-initio perspective. 
We study the Mi phase of VO2 using first-principles cal- 
culations as a function of static disorder with a state of 
the art linear scaling DFT method [15]. The capabil- 
ity of linear scaling DFT to describe large super-cells, 
containing several hundreds of atoms, is necessary to 
comprehensively tackle the issue of disorder. We extend 
our DFT calculations with the dynamical mean-field the- 
ory (DMFT) approximation [16, 17] in order to refine 
the description of the strong correlations induced by the 
3d subshell of the vanadium sites (for more details see 
the supplementary material). Throughout this work we 
used typical values for the screened Coulomb interaction 
(U = 4 eV) and Hunds coupling (J = 0.68 cV) [18, 19], 
and our calculations were carried out for 324 atom super- 
cells (108 V atoms) and 768 atom super-cells (256 V 
atoms) at fixed temperature T = 189 K. All orbitals 
are defined in the local coordinate system [20] associated 
with the vanadium atoms. 

We first discuss single-site DMFT calculations for 
paramagnetic VO2 . The dependence of the spectral func- 
tion on the on-site repulsion Ud is shown in Fig. 1. We 
find that the Mi phase of VO2 is metallic for Ud = 4 eV 
and that there is a large spectral weight present at the 
Fermi level. Hence, V0 2 is described by DFT+DMFT 
as a charge transfer correlated paramagnetic metal, with 
a moderate mass renormalization m* /m — 1.35, of the 
same order as the mass renormalization in the rutile 
phase obtained by other groups (m*/m = 1.8 from 
Ref. [3] and m* /m = 1.51 from Ref. [10]). The large spec- 
tral weight at the Fermi level in Fig. 1 is of predominantly 
d xy character, the contribution from the e g orbitals being 
negligible: the spin-independent orbital densities at the 
Fermi level p a {e F ) are 0.02,0.02,0.19,0.25,0.28 for, respec- 
tively, d x 2_ v 2,d 3z 2_ r 2,d yz ,d xz ,dxy symmetry, which indi- 
cates a strong selection of the ta 9 orbitals at the Fermi 
level in agreement with the orbital selection scenario ar- 
gued long ago by Goodenough [7]. Notably, we find that 
the dynamical correlations, described by the imaginary 
part of the self-energy, also suggest that the d xy orbital 
is the most correlated orbital, whereas the e g states are 
weakly correlated (for more details see Fig. 2 and Fig. 3 
of the supplementary material) . 

However, our results deviate from the description of 
Goodenough. In particular, we obtain a vanadium 3d 
sub-shell filling of n = 3.15 electrons from DFT, much 
larger than the 3d 1 configuration of the ionic picture. 
We emphasize that we used a set of local Wannier or- 
bitals, variationally optimized during the energy mini- 
mization carried out in the DFT calculations [21], which 
renders the calculation of the electronic density very re- 
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FIG. 1: ZSA boundary: (Color online) Dependence of the 
spectral function p(ui) with respect to the Coulomb repulsion 
U d . 



liable. Larger 3d orbital occupations in VO2 than the 
single electron have been reported in earlier DFT calcu- 
lations (LSDA+t/ finds n = 2.48 e [6]). We note, fur- 
thermore, that similar occupancies are obtained for the 
R phase, both by experimental measurement (n = 1.78 e 
from XAS [6]) and DFT calculations (LSDA+C/ finds 
n = 2.31 e [6] and LMTO-ASA gives n = 3.35 e [22]). 

For larger Ud, we find that the spectral weight at 
the Fermi level shrinks, and we obtain an insulator for 
Ud = 25 eV, placing VO2 well below the ZSA bound- 
ary U%, estimated between 21 eV and 25 eV. Since we 
obtain a paramagnetic insulator only for Coulomb re- 
pulsions of the order of the atomic limit, this strongly 
suggests that Mott correlations alone are not responsible 
for the insulating state of the M x phase, as also suggested 
by early LDA calculations [20, 23] which failed to repro- 
duce the insulating state. We conclude that a large 3d-3d 
Coulomb interaction alone is not sufficient to generate a 
large band gap for VO2. Finally, we also explored the 
dependence on the Hund's coupling J and found no sig- 
nificant change in the mass renormalization by varying 
J between 0.3 eV and 1.2 eV, for fixed Ud = 4: eV, al- 
though increasing J enhances the mass renormalization 
for U d = 8 eV and U d = 16 eV. 

If Coulomb correlations alone cannot lead to insulat- 
ing behavior, perhaps the inevitable disorder due to im- 
perfections of the crystal, or self-trapping due to strong 
electron-phonon coupling could be relevant. Hence we 
applied a random three-dimensional Gaussian displace- 
ment to both the V and O atomic sites. The Gaussian 
width 5 characterizes the amplitude of the disorder. We 
also note that a large number of uncontrolled and uncor- 
rected experimental parameters, such as the variation of 
the atomic positions due to the shear induced by dislo- 
cations, are described at the theoretical level by a Gaus- 
sian distribution (central limit theorem). Notwithstand- 
ing, we use a single representation of disorder for each 
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FIG. 2: Anderson and Mott localization: (Color 
online) a) Spectral function p of paramagnetic VO2 in the 
presence of Gaussian disorder S. b) Averaged quasi-particle 
weight Zd with respect to the repulsion Ud for zero (5 = A) 
to large disorder (S = 0.5 A). Calculations including a single 
O vacancy in the S = A case are also shown for compari- 
son (open squares), c) Distribution of the local quasi-particle 
weight Zd,i for S = 0.1 A. d) Isosurface of the real space rep- 
resentation of the Fermi density for disorder 5 — 0.3 A. The 
large (small) sphere denotes V (O) atoms along the rutile axis. 



amplitude 6, Gaussian disorder satisfying the ergodicity 
theorem such that the local average over different config- 
urations is equivalent to the average over space for large 
enough super-cells. 

The spectral function for disordered VO2 is shown in 
Fig. 2. a. Although the spectral weight at the Fermi level 
is suppressed with increasing disorder (reflecting charge 
localisation) the system remains metallic up to the largest 
physical amplitudes of disorder. The effect of the local- 
ization induced by disorder is also observed in the aver- 
aged quasi-particle weight Zd (Fig. 2.b) and in the spa- 
tial distribution of the local quasi-particle weight Z^i 
(Fig. 2.c), which clearly shows that the disorder gener- 
ates regions in the crystal with strong localization, which 
coexist with metallic parts of the crystal where the local- 
ization has only a weak effect. These droplets of strongly 
correlated Fermi liquid generate a larger mass renormal- 
ization m* /m on average, as observed in the decrease of 



the averaged quasi-particle weight (Zd = m/m*) as the 
disorder increases (Fig. 2.b). 

The localization effect can be understood in a simple 
picture: when the O atoms move closer to the V atomic 
site, the static charge repulsion induces a larger charge 
transfer energy A = — e p , which enhances the strength 
of the correlation locally (the repulsion U of the one band 
Hubbard model translates into the charge transfer energy 
in d-p theories [13]). This effect is illustrated in Fig. 2.d, 
where we show an isosurface of the real-space representa- 
tion of the Fermi density p(ep, r) for one of the V chains 
along the rutile axis for 5 = 0.3 A, where large (small) 
spheres denote V (O) atomic sites. The V atom high- 
lighted by a star has two very near oxygen neighbors, 
which is expected to induce a larger charge transfer en- 
ergy. The latter results in a transfer of charge from the 
vanadium site to one of its oxygen neighbors (indicated 
by an arrow). The subtle interplay between the local- 
ization induced by the disorder (Anderson-like) and the 
localization induced by strong correlations (Mott-like) is 
captured by the DFT+DMFT methodology. 

We now move to the non-local cluster cellular 
DMFT calculations (c-DMFT). The non-local correla- 
tion present in cluster DMFT drives VO2 to an insulator 
(Fig. 3. a), in agreement with earlier DMFT calculations 
using model Hamiltonians [10] and we obtain a gap of 
~ 0.6 eV in agreement with both the latter and the ex- 
perimental value [9]. We did not observe any finite size 
effect, and the Peierls gaps obtained extracted from the 
324 and 768 atom super-cells are identical. We find that 
large disorder quenches the Peierls state for S > 0.1 A. 
Very interestingly, the insulating Peierls state survives 
for moderate disorder 6 — 0.1 A, although the gap is 
reduced down to ~ 0.3 eV. We also carried out cluster 
DMFT calculations for the case of a single O vacancy: 
the O vacancy creates a mid-gap state (inset of Fig. 3. a), 
spatially localized in the center of the three V atoms sur- 
rounding the vacancy, as illustrated by the real space 
representation of the Fermi level density (Fig. 3.b) but 
does not strongly affect the band edges. In conclusion, 
our results suggest that the Peierls instability in V0 2 is 
very robust, surviving external perturbations such as re- 
duction of the long-range crystallographic order or local 
impurities. 

The imaginary part of the self-energy of the dynamical 
Peierls singlet is shown in Fig. 3.c. We observe that the 
gap is mainly induced by dynamical correlations in the 
d xy orbital, which exhibit a pole at the Fermi level. In 
our view, the dynamical V-V dimers generate a Mott 
instability (the mechanism may be thus termed Peierls- 
Mott). In particular, the spectral weight (inset) shows 
that the cluster DMFT almost entirely depletes the d yz 
orbital, leaving two electrons equally shared between the 
d xz and d xy orbitals. The lobe of the latter orbitals point 
towards the rutile axis, whereas the d yz orbital is oriented 
perpendicular to this direction and thus the latter does 
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FIG. 3: Mott-Peierls instability: (Color online) a) Spec- 
tral function p obtained using cellular cluster DMFT (c- 
DMFT) calculations without (S — A) disorder for moderate 
(324 atom) and large (768 atom) super-cells. Calculations 
for disordered VO2 (5 — 0.1 A) and for a single O vacancy 
are also shown for comparison. Inset: Enlargement of the 
low energy scale, the gap of ~ 0.6 eV is shown. The va- 
cancy introduces a mid-gap state, highlighted by the arrow, 
b) Isosurface of the real space representation of the charge 
density at the Fermi level for calculations for an O vacancy 
(star). The large (small) sphere denotes V (O) atoms along 
the rutile axis (horizontal direction), c) Imaginary part of the 
self-energy and (inset) imaginary part of the Green's function 
for 5 = A. 



not contribute strongly to the orbital bonding within a 
V-V dimer. 

Interestingly, therefore in our picture we find that 
two electrons per V atom lie in bonding orbitals, lead- 
ing to a strong Mott dynamical divergence in the self- 
energy (Peierls-Mott). This contrasts with the picture of 
Ref. [10], where a single electron on each V is of bonding 
character and the repulsion Ud drives the bonding or- 
bitals to a singlet configuration, following the early pro- 
posal of Sommcrs and Doniach [24]. In the latter picture, 
the repulsion energy may be dramatically reduced by the 
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FIG. 4: Optical spectra: (Color online) a) Theoretical op- 
tical conductivity along the rutile axis (solid lines) and along 
the perpendicular direction (dashed lines) calculated with cel- 
lular DMFT. Experimental data for poly crystalline films [26] 
(squares) and thin films (diamonds) [27] are shown for com- 
parison, b) Optical conductivity along the rutile axis obtained 
by cellular DMFT for disordered VO2 for various disorder am- 
plitudes S. 



formation of the singlet state, manifested in the fact that 
the low-frequency behavior of the on-site component of 
the self-energy, that associated with the ai g orbital, is lin- 
ear in frequency, as opposed to a Mott insulator in which 
S" diverges. In our picture, orbital rehybridization acts 
to deplete the d yz orbital, leaving 2 e in two orbitals per 
V site, in turn generating a Mott instability which cre- 
ates a pole in the local self-energy (Peierls-Mott). We 
note that this result reconciles the Peierls scenario with 
recent XAS [6] and photoemission spectroscopy measure- 
ments [25], which hint at an occupation of n ss 2 for the 
3d sub-shell in the Mi phase. 

Finally, the optical conductivity calculated using clus- 
ter DMFT (Fig. 4. a) is in qualitative agreement with 
experimental data obtained for polycrystalline films [26] 
and thin films [27]. We note that the optical gap is not 
dramatically affected by a moderate degree of disorder 
5 = 0.1 A (Fig. 4.b). For large disorder, however, VO2 
is a bad metal and we note, in particular, that no Drudc 
peak is obtained in the optical conductivity, and that 
the disorder induces strong oscillations in the optical re- 
sponse for the infrared frequency range 10 < 1 eV. 

In conclusion, we have carried out linear scaling 
first principle calculations, in combination with clus- 
ter DMFT, on VO2, both with and without disorder. 
We find that VO2 is a Peierls-Mott insulator, and that 
the ZSA boundary of the paramagnetic insulator is ob- 
tained only for unrealistic values of the Coulomb repul- 
sion (Ud ~ 25 eV). We propose a new mechanism for the 
insulating Mi phase of VO2 based on an orbital selec- 
tive Mott transition, assisted by the Peierls distortion: 
the Peierls instability involves an orbital selection, and 
bonds the d xy and d xz orbitals along the rutile axis, fill- 
ing each orbital with one electron, and in turn gener- 
ates a Mott instability. This scenario may be described 
as Peierls assisted orbital selective Mott transition. Fi- 
nally, we demonstrated that the Peierls phase survives 
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moderate Gaussian disorder (5 — 0.1 A), and hence our 
picture accounts for the observation of the MIT in the 
experimentally realistic, disordered system [5]. Whether 
the critical temperature associated to the Mott instabil- 
ity is lower than the structural phase transition temper- 
ature (SPT), and whether this picture accounts for the 
recent observation of metallic responses in VO2 in the Mi 
phase near the SPT [28] remains to be seen. Finally, we 
found that oxygen vacancies induce a localized mid-gap 
state, leaving the band edges unaffected, shedding some 
light on thin-film measurements where substrate strain 
can induce stoichiometric modification [8]. Our results, 
combining lattice disorder and a powerful method for de- 
scribing non-local, dynamical correlation, open up new 
frontiers for first principles materials design under under 
realistic experimental conditions. 
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METHODOLOGY 

Kohn-Sham density- functional theory (DFT) [1, 2] has 
been shown to be capable of making accurate predictions 
for many materials. DFT combines high accuracy and 
moderate computational cost, yet the computational ef- 
fort required to perform calculations with conventional 
DFT approaches is still non negligible: it increases with 
the cube of the number of atoms, rendering them unable 
to routinely tackle problems comprising more than a few 
hundred atoms, even on modern large supercomputers. 
Since the pioneering work of the Nobel laureate Walter 
Kohn [3] , it has been known that it is possible to reformu- 
late DFT to exploit the short-ranged nature of quantum 
mechanics, so that the computational cost scales linearly 
which would, in principle, allow calculations with many 
thousands of atoms. The ONETEP approach [4-6] , used 
in this work, is an example of the practical realization of 
linear-scaling DFT and DFT+J7, and is particularly no- 
table for its accuracy equivalent to a plane- wave method, 
by means of its in situ optimization of a set of local Wan- 
nier orbitals with a respect to a systematically improv- 
able basis. 

The linear scaling performance of ONETEP allowed 
us to study large atomic super-cells. Specifically, we in- 
vestigated both a moderate supercell of 324 atoms for 
VO2 (108 V atoms) and a large supercell of 768 atoms 
(256 V atoms). We used, for the crystallographic struc- 
ture of the Mi phase of VO2, the lattice constants and 
atomic positions obtained for powder VO2 (a = 5.743 A, 
b = 4.517 A, c = 5.375 A, and (3 = 122.6°) since the 
lattice parameters of thin films vary with charge local- 
ization and substrate- induced built-in strain [7]. The 
electronic core states were represented by GGA scalar- 
rclativistic pseudo-potentials generated with the OPIUM 
package [8]. The (3s, 3p, 3d, 4s) and (2s, 2p) atomic shells 
were respectively treated as valence states for the V and 
O atoms. Hence, we used 10 orbitals to describe each 
V atom and 4 to describe each O atom, giving a total 
of 1944 orbitals for the 324 atom super-cell and 4608 for 
the 768 atom super-cell. 



ONETEP [4] is at the cutting edge of developments in 
first principles calculations. However, while the funda- 
mental difficulties of performing accurate first-principles 
calculations with linear-scaling cost have been solved, one 
remaining problem of DFT approaches in general is that 
conventional approximations to the exchange-correlation 
(XC) functions fail in describing many compounds where 
strong correlations are present. In such cases, DFT of- 
ten predicts moments and energetics that are qualita- 
tively inconsistent with experiment, or fails to describe 
the insulating state and converges to a metallic state in- 
stead. DFT+DMFT is a computationally moderately 
expensive method for refining the description of on-site 
Coulomb interactions provided by such XC functional 
and, hence, for extending the range of applicability of 
DFT to strongly-correlated materials. 

The principle of DFT+DMFT is to separate the dom- 
inant subspace of the system, which is well-described by 
conventional XC functionals due to its delocalized and 
free-electron character, from a subset of localized sub- 
spaces, so-called correlated sites, which are strongly cor- 
related and not adequately described by the XC func- 
tionals. The XC correlated sites are then explicitly aug- 
mented with screened Coulomb and Hund's coupling in- 
teractions within these sites, together with a double- 
counting term to correct for the component already in- 
cluded within the DFT XC functional. 

The correlated sites are typically spanned by a set of 
3d and/or 4f atomic-like orbitals, which define a set of lo- 
calized Hubbard projectors [9]. Solutions of appropriate 
orbital symmetry of the hydrogenic Schrbdinger equa- 
tion, such as atomic-like or linear muffin-tin orbitals, are 
the most common choice [10-12] for Hubbard projectors. 
The vanadium-centered correlated sites were delineated 
in our calculations using hydrogenic 3d orbitals as Hub- 
bard projectors, characterized by the canonical Clementi- 
Raimondi [13] effective charge of Z ~ 8.9829 for this 
atomic sub-shell and species, determining their spatial 
diffuseness. 

In the implementation of the linear-scaling DFT 
method used in this work, we work with the single- 
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particle density-matrix, which is expressed in a separable 
form [14, 15] [45]: 

p(r,r')=^0 Q (r)K"%(r'), (1) 

a/3 

where {<j) a } is a set of localized nonorthogonal gener- 
alized Wannicr functions (NGWFs) variationally opti- 
mized during the energy minimization carried out as 
part of the DFT calculations [16]. The density kernel 
K a/3 = {(j) a \p\<j)P) is the representation of the single- 
particle density operator p in terms of the contravari- 
ant NGWF duals {<fi a }, defined as those which satisfy 
the relationship (cj) a \(f)p) = 8 a p. The NGWFs are, in 
turn, expanded in terms of a systematic basis of Fourier- 
Lagrange, or psinc functions [17]. The size of this basis 
is determined by a plane-wave kinetic energy cutoff (a 
minimum of 822 eV for all reciprocal lattice vectors in 
our calculations). The DFT energy functional is itera- 
tively minimized with respect to both the density kernel 
and the NGWF expansion coefficients using a minimiza- 
tion scheme which is detailed in Ref. [18]. Therefore, the 
NGWFs are a readily accessible set of localized orbitals 
(spatially truncated to atom-centered spheres of radius 
6.6 A) which are calculated with linear-scaling compu- 
tational cost. No additional approximation such as the 
density kernel truncation was used in this work, and the 
energy convergence was tested against the plane-wave ki- 
netic energy cutoff and the spatial truncation of the NG- 
WFs. 

Once the fully converged energy minimization of VO2 
was carried out, the full Green's function was computed 
in the finite temperature Matsubara representation [46] 
from the readily available DFT Hamiltonian H: 

G afs (tw n ) - ((iw n + n)S aP - H afj - E^r 1 , (2) 

where /1 is the chemical potential (fixed at the mid-energy 
between the last occupied state of the ONETEP calcu- 
lation and the first empty state), S is the full overlap 
matrix between the NGWFs so that <j) a — S a p4>^ an d 
S is the self-energy matrix computed by the DMFT al- 
gorithm. We note that the computation of the Green's 
function requires large matrix inversion (in particular, for 
VO2, the matrix to invert is of size 1944 x 1944 for the 
324 atom super-cell and 4608 x 4608 for the 768 atom 
super-cell. These super-cells allowed us to describe spa- 
tially dependent perturbations, such as static disorder 
(see Fig. 1), and also to consider oxygen vacancies [47]. 
Matrix inversion was carried out for 100 Matsubara fre- 
quencies to provide converged sampling of the Green's 
function. To achieve reasonable computational time, we 
performed the matrix inversion on recently-developed 
graphical computational units (GPUs) using a home- 
made parallel implementation of the Cholesky LU decom- 
position with the CUDA programming language. The 
calculations were carried out using 8 GPUs in parallel, 




FIG. 1: Super-cell for DFT+DMFT calculations: 

(Color online) Super-cell used for the 324 atom DFT+DMFT 
calculations for disordered VO2 with S = 0.1 A. 



achieving a total of 4000 GFLOPS (giga floating-point 
operations per second), allowing a tremendous speed-up 
of Green's function calculation (the computational time 
needed to carry out the matrix products and matrix in- 
versions involved in the calculation for a 100 Matsubara 
frequencies was of the order of a few minutes for the 
largest super-cell). The matrix multiplications involved 
throughout the DMFT calculations were also performed 
using the CUDA architecture. 

The projected Green's function G of the vanadium 3d 
correlated sites is obtained as 

W = W<&G a ' ) (io,„) t$2„ (3) 

where I runs over the Vanadium atoms, m and m' over 
the five atomic 3d orbitals used as orthonormal Hub- 
bard projectors (in real cubic-harmonic notation: d x 2_ y 2 , 
d 3z 2_ r 2, d yz , d xz , d xy ), a and /? are the indices for 
the NGWFs, and the matrices Vam — ((f> a \<Pm) and 
Wml — (ipm\<f>a) are the overlap between the NGWFs 
{4> a } and the hydrogenic orbital {ipi^} at the Vanadium 
site I and with orbital index m. Note that the orthonor- 
mality of the Hubbard projectors implies that we do not 
need to distinguish them from their duals (the generaliza- 
tion to nonorthogonal projectors is discussed in Ref [19]). 

The strong correlation acting within each Vanadium 
site is described by the conventional Slater-Kanamori 
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form of the interaction term [20, 21], namely 



Hint — Ud 



m>m / 



- J ( 2S ™Sm> + (<4t44 d m'T d m'i)), 

m>m / 

for each correlated site and with no inter-site interactions, 
where the first term describes the effect of intra-orbital 
Coulomb repulsion Ud and the second term describes the 
inter-orbital repulsion U' d . The latter is renormalized by J 
in order to ensure rotational invariance of the interaction 
(for a review see Ref. [22]). The third term is the Hund's 
rule coupling, described by a spin exchange of amplitude 
J. S m denotes the spin of orbital 771, S m — ^d^gCTgs'dms' , 
where <r is the vector of Pauli matrices indexed by s and 
s'. Generally, the hierarchy of interactions is Ud > U' d > 
J. Ud and J were set to typical values for the screened 
interactions for VO2, namely Ud = 4c V and J = 0.68eV 
[23, 24]. 

The interaction term (4), combined with (2) and (3), 
defines a local problem at each Vanadium site I of the 
lattice, which we solve with the DMFT algorithm [25]. In 
the single site DMFT approach, the projected self-energy 
matrix £( J )( J ) is assumed to be local, i.e., 



fj(0(J) = 



(5) 



for non-local Vanadium site indices (I 7^ J), unlike the 
cluster cellular DMFT [26] (c-DMFT) approach where 
spatial non-local contributions from vanadium sites be- 
longing to the same dimer are included. 

The self-energy is obtained by solving an auxiliary An- 
derson impurity problem (AIM) for each V atom (V-V 
dimer) in the single DMFT (cluster DMFT). The hy- 
bridization matrix A in the orbital subspace (d x 2_ y 2, 
d 3z 2_ r 2, d yz , d xz , d xy ) of the AIM problem is defined, 
for each site, as: 

A(^„) = (7w„+ M )6-S-E im P-G- 1 , (6) 

O is the local projected overlap matrix, which is related 
to the non-orthogonality of the DFT basis set (NGWFs) 
[27-29] and O = 1 in the particular case of orthonormal 
Hubbard projectors in the form of Kohn-Sham Wannicr 
functions. In particular, we find that O is related to the 
NGWF's overlap matrix by 



(7) 



The static part of the hybridisation function W mp = 
A(iw„ — > 00) defines the energy levels of the im- 
purity in the absence of hybridisation (atomic limit). 
Both O and E lmp define the impurity problem and 
they are constructed such that the hybridization ma- 
trix A(i(x>„) exhibits the correct physical decay cx 



l/iuj n at large frequency (we ensured that 0^ m , — 
lim^^co G 1 (iuj) juj to a high precision, as required 

by this condition). The self-energy £ is obtained by 
solving the Anderson impurity problem defined by the 
hybridization (6) and the interaction Hamiltonian (4). 

The AIM is solved using two different solvers: (1) a fi- 
nite temperature Lanczos algorithm [30-32] , and some of 
the results were cross-checked and benchmarked with (2) 
a continuous-time Monte Carlo solver (CTQMC) [33, 34]. 
The two solvers both suffer from approximations but 
complement each other: the Lanczos solver uses a fi- 
nite discretization of the hybridization (6) and suffers 
from finite size effects. The CTQMC solver suffers from 
the so-called fermionic quantum sign problem when the 
off-diagonal elements are considered, hence in this work 
we neglect within the CTQMC implementation the off- 
diagonal elements of the hybridization matrix (6). We 
note, however, that the amplitude of the off-diagonal el- 
ements depends strongly on the choice of the localized 
Hubbard projectors, and, in particular, they depend on 
the choice for the local orientation of the local (e^ey^) 
axis. The off-diagonal elements of the hybridization and 
Green's function prove to be small for the vanadium 3d 
set employed in this work when the local coordinates de- 
fined for VO2 are used (see discussion hereafter). 

Once the local self-energy £ is obtained, it is up-folded 
to the space spanned by the Kohn-Sham orbitals (equiv- 
alently the space spanned by NGWFs), in the separable 
form 



/ j am^mm' rr m'f3 
I, J 



(8) 



As pointed out recently in [35], the separable form of 
the projectors of the self-energy enforces the causality 
of the up-folded self-energy (the local self-energy in the 
Hubbard projector basis is causal by construction, so 

that Tr m 'm {iuj n ) < 0, V(iw n ,m, /), since it is obtained 
by solving the AIM in the local orthogonal hydrogenic 
basis). We note that the causality is a necessary re- 
quirement to obtain physical observables. We carefully 
checked that all of the obtained self-energies were causal, 
since numerical error might still induce some unexpected 
fluctuation in the imaginary part of the self-energy). 

One particular difficulty arises when one wants to con- 
struct a non-interacting Hamiltonian H from the DFT- 
LDA ground-state density. Since the LDA already con- 
tains the influence of the Coulomb interaction to a certain 
degree, the problem of double counting of these contribu- 
tions by the interaction Hamiltonian (4) arises. In order 
to avoid this double counting, one has to subtract off the 
interaction contributions from the LDA. Unfortunately, 
the precise form for a particular set of orbitals is not 
known and the best one can do is to account for these 
contributions in an averaged way. In this work, we used 
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the canonical form of the correction E dc , given by 

E dc = Uf (n d - 0.5) -i( na -l), (9) 

where rid is the total occupancy of the vanadium 3d sub- 
space, as defined by the Hubbard projectors, and U% v 
is the averaged repulsion related to the intra-orbital and 
inter-orbital repulsion, computed as follows [36]: 



US 



Ud + 2 (iV dcg - 1) U' d 
2N dcg - 1 



(10) 



and N deg is given by the number of orbitals. 

The upfoldcd self-energy, in the covariant NGWFs rep- 
resentation, is hence corrected by E dci giving 

i.j 

The double-counting correction avoids spurious artifacts, 
such as a shift of the subspace occupancies when the 
Coulomb interaction Ud is increased from zero. In this 
work, we did not find any significant deviation from the 
DFT density in the DFT+DMFT scheme, providing a 
strong indication that the double-counting term (9) cor- 
rectly accounts for interactions already included in the 
DFT-LDA Hamiltonian H. 

For the particular case of V0 2 in the Mi phase, due 
to the buckling of the V atoms, it is convenient to carry 
out the DMFT in a rotated basis, making use of a local 
coordinate system associated with the vanadium sites. In 
particular, the rotation used in this work is such that the 
local e x axis is directed towards the rutile axis and the 
local e y axis points towards the in-plane oxygen. The 
vanadium- apical oxygen vector forms a 90° angle with 
the rutile axis [37]. This allows us to obtain an hybridiza- 
tion matrix (6), which is as diagonal as possible. This is 
crucial for the CTQMC algorithm, since it neglects the 
off-diagonal terms. The rotation is carried out in orbital 
space [38] and translates for the projected Green's func- 
tion as follows: 



Qrot = utGU, 



(12) 



where U is the 5x5 rotation matrix in cubic harmonic 
space, corresponding to a real space rotation R( 3 \ and 
G rot is the rotated Green's function. The same holds for 
the projected self-energy obtained by solving the AIM 
model S rot , which is rotated back to the original system 
of coordinates: 



s = us rot u t . 



(13) 



Equations (2), (3), (6) and (11) are iterated until DMFT 
convergence is obtained. Real frequency quantities are 
readily available from the Lanczos solver and provides us 



access to the physical observables of the system. In par- 
ticular, the density of states is available from the retarded 
full Green's function: 



(14) 



where the NGWF-resolved spectral density p a p (lj) is de- 
fined by the retarded Green's function obtained at the 
DMFT self-consistency Gr as: 



n a/3 



2in 



fa/3 , 



(15) 



A very common experimental probe is optical spec- 
troscopy, which is accessed in this work through the 
optical conductivity. The optical conductivity in 
DFT+DMFT is computed by using the Kubo formula of 
linear response theory in the no-vertex-corrections [48] 
approximation [39] , via the expression 



f(u/-u)-f{«/) 



2ire 2 h 

n 

x (p afj (uj'-w)v 01 p^ 5 (u')v Sa ), 



(16) 



where the factor two accounts for the spin degeneracy, 
is the simulation-cell volume, e the electron charge, 
h the reduced Planck constant, / (u>) is the Fermi-Dirac 
distribution, and the optical conductivity is obtained by 
the convolution of charge density matrix operator p a P , 
defined in equation (15), and the bare vertex v a( g is given 
for a non-orthogonal basis by [40] 



w a p = p/m e = -ih/m e (<j) a \V\<pp). 



(17) 



We emphasize that the so-called Peierls substitution [39] 
approximation was not used here. 

Another interesting physical quantity probed in this 
work is the quasi-particle weight of a local vanadium 3d 
subspace Z™^ where the index i runs over the Vanadium 
site and m over the d orbitals. We used the discrete 
Eliashberg estimate to obtain the quasi-particle weight: 



v/(i)(») 



ZT, = 1 - 



(18) 



where u>i is the first Matsubara frequency. The quasi- 
particle weight Z d ,i of the Vanadium site i is estimated 
by the Fermi liquid specific heat of each independent d 
orbital, that is [41]: 



li oc 



L dA 



where ep is the Fermi energy, and it follows: 



(19) 



(20) 
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The effect of the localization induced by disorder can be 
quantized by the average Zd, which is obtained along the 
same lines: 

Z«i = £>M/I>v (21) 

i i 

Zd is related to the mass renormalization {Z& — m/m*). 



COMPARISON BETWEEN THE ED AND 
CTQMC IMPURITY SOLVERS 

In this section, we extend the discussion of the single 
site DMFT results for V0 2 . We first discuss the Mat- 
subara representation of the Green's function (the real 
and imaginary parts of the Green's function are shown 
respectively in Fig. 2. a and Fig. 2.b.) We obtained, in 
particular, that the orbital density at the Fermi level 
n = -G" (icj n = 0) /tt obtained by CTQMC (open sym- 
bols) is largest for the d xy orbital. Interestingly, we find 
that the e g states are almost empty and do not have a 
significant weight at the Fermi level. The case for such an 
orbital selection process was argued long ago by Goode- 
nough [42] : VO2 can be regarded as a 3d 1 system in the 
simple ionic picture (with full charge transfer from the 
V 3d 3 4s 2 to the two O 2p 4 configurations). The low en- 
ergy states near the Fermi level are of strong vanadium 
3d character and each vanadium atom is surrounded by 
an oxygen octahedron, such that the crystal field splits 
the 3d states into t2 9 and e g bands. Since the structure 
is not cubic, however, the t 2 g states are further split into 
e^g and &i g states. In this simple ionic picture, the re- 
maining 3d electron partially occupies the ai g band and 
the system remains metallic in the rutile phase. In the 
Mi phase, pairs of vanadium atoms along the rutile axis 
form dimers. This Peierls distortion causes strong hy- 
bridization between the ai ff orbitals of the two vanadium 
atoms, allowing the bonding state to fully fill and open- 
ing a gap between the bonding ai g and the unoccupied 
7r bands. 

While the t^g orbitals point between the oxygen sites, 
the d 2; 2„j / 2 and ds z 2- r 2 orbitals point towards the oxy- 
gen sites, giving strong cr-bonds with the O 2p orbitals: 
because of the high electronegativity of oxygen (3.5) 
compared to vanadium (1.6), the bonding d 32 2_ r 2 and 
d x 2_ y % orbitals are of mainly O 2p character and hence 
the electrons are more localized to the O 2p than the V 3d 
orbitals; an almost ideal ionic bond. This suggests that 



the most localized orbitals in the 3d shell are the t2 S , 
and that polarized X-ray spectroscopy mainly resolves 
t 2g orbitals at the vanadium sites. 

This scenario is also supported by the CTQMC data 
(Fig. 2.a,b), which shows that dynamic correlations ef- 
fects are not important for the e 9 orbitals. We obtained, 
from DFT+DMFT solved with the CTQMC, total densi- 
ties of (0.68,0.73,0.52,0.52,0.68) respectively in the com- 
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FIG. 2: Matsubara representation of the Green's 
function: (Color online) Real (a) and imaginary part (b) 
of the Green's function G in the Matsubara representation 
for paramagnetic VO2. The calculations were performed us- 
ing the CTQMC impurity solver (open symbols) including 
all five Vanadium d orbitals (e 9 and t2 g ) are compared with 
the calculations performed with the ED impurity solver (filled 
symbols), which includes only the three ti B orbitals. 



ponents of the 3d subshell (d x 2_ y 2 ,d3 Z 2_ r 2 ,d yz ,d xz ,d xy ). 
Thus, there is n = 1.4 e in the e g and 1.7 e in the t2 S 
orbitals, in broad agreement with the n = 1.83 e occu- 
pation of the 3d shell extracted from XAS measurements 
[43] and n = 1.9 e obtained using photoemission spec- 
troscopy [44]. 

The real and imaginary parts of the self-energy are 
shown in Fig. 3. a and Fig. 3.b. The dynamical correla- 
tions, described by the imaginary part of the self-energy, 
show that the d xy orbital is also the most correlated or- 
bital, whereas the e g states are weakly correlated. The 
latter is also supported by the weak frequency depen- 
dence of the e ff states in the real part of the self en- 
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FIG. 3: Matsubara representation of the self-energy: 

(Color online) Real (a) and imaginary part (b) of the self- 
energy £ in the Matsubara representation for paramagnetic 
VO2 • The calculations performed with the CTQMC impurity 
solver (open symbols) including all five Vanadium d orbitals 
(e s and £29) are compared with the calculations performed 
using the ED impurity solver (filled symbols), which includes 
only the three t^g orbitals. 
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FIG. 4: Real space representation of the scattering 
rate for cellular DMFT calculations: (Color online) 
Contour plot of the real-space representation of the scattering 
rate at the Fermi level — E" (r, r) obtained by cellular DMFT 
(c-DMFT) shown for a [010] cut of the supercell-cell with 
zero (a) and finite disorder (5 = 0.1 A) (b). c) Real-space 
iso-surfaces of the scattering rate at the Fermi level for the 
case of the O vacancy (star). The large (small) spheres denote 
V (O) atoms. The rutile axis is oriented along the horizontal 
direction in all the calculations above. 



(filled symbols in Fig. 2.a,b and Fig. 3.a,b) and CTQMC 
(open symbols in Fig. 2.a,b and Fig. 3.a,b) is remark- 
able for both the Green's function (Fig. 2.a,b) and self- 
energy (Fig. 3.a,b). Finally, we obtain also a very good 
agreement for the quasi-particle weight Zd obtained with 
the ED (Z d = 0.77) and CTQMC (Z d = 0.73) solvers. 
The data shown in the paper and supplementary mate- 
rial were obtained with the ED solver when not specified 
otherwise. 



REAL SPACE REPRESENTATION OF THE 
SELF-ENERGY OBTAINED BY C-DMFT 



ergy (Fig. 3. a). Since the e g orbitals do not exhibit any 
strong dynamical correlation, it is a good approximation 
to treat them at the Hartree level and thus we applied 
the ED solver only to the t2 g orbitals. We tested this 
approximation by comparing the densities obtained with 
the two methods. In particular, we obtain by using this 
approximation and the ED solver a total density in the 
t2 g orbital n = 1.69 e, in very good agreement with the 
CTQMC calculation, which finds a density of n = 1.7 c 
in the t2 9 orbitals, and includes all the five d orbitals 
within the calculation. The comparison between the ED 



In this section, we extend the discussion regarding 
the cellular c-DMFT calculations under the presence of 
external perturbations, such as disorder and vacancies. 
Contour plots of the real-space representation of the lo- 
cal part of the scattering rate, — S"(r, r), obtained by 
cellular DMFT (c-DMFT), for a [010] cut of the super- 
cell, are shown for zero and finite disorder in Fig. 4. a 
and Fig. 4.b, respectively. Interestingly, we find that 
a moderate disorder of S = 0.1 A enhances the dimer 
bonding character by reducing the V-V distance among 
some of the dimers and thus increasing the self-energy. 
This is illustrated by the real-space representation of the 
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scattering rate, in the [010] plane, ol the super-cell for 
5 = A (Fig. 4.a) and S = 0.1 A (Fig. 4.b). Signifi- 
cantly, although the bonding character of some dimers 
is inevitably diminished, the integrated spectral function 
(Fig. 3. a in the paper) shows that the Mott localization 
effect is strong enough to sustain an insulator, such that 
no spectral weight is present at the Fermi level. 

A very different effect is observed for the case of a sin- 
gle O vacancy. At 5 = A, this defect does not affect the 
band edges but instead induces a mid-gap state with a fi- 
nite weight at the Fermi level. The real-space iso-surfaccs 
of the scattering rate show (Fig. 4.c) that the vacancy 
breaks a dimer (dimer C) and enhances the bonding of 
the nearest neighbor dimer (D). In particular, the strong 
mass renormalization associated with the Mott instabil- 
ity prevents the charge to flow freely between neighbour 
dimers, such that the mid-gap state remains essentially 
local in space. We argue that this mid-gap, localized 
state is relevant to the thin-film deposition of V0 2 lay- 
ers, where it was recently found that the strain of the 
substrate might induce subtle deviation from the canoni- 
cal stoichiometry, and moreover that the metal-insulator 
transition was strongly varied over different samples [7]. 
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